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Large scale molecular dynamics simulations are performed to study the steady state yielding 
dynamics of a well established simple glass. In contrast to the supercooled state, where the shear 
stress, o, tends to zero at vanishing shear rate, 7, a stress plateau forms in the glass which extends 
over about two decades in shear rate. This strongly suggests the existence of a finite dynamic 
yield stress in the glass, a + (T) = a(T;7 — > 0) > 0. Furthermore, the temperature dependence of 
a + suggests a yield stress discontinuity at the glass transition in agreement with recent theoretical 
predictions. We scrutinize and support this observation by testing explicitly for the assumptions 
(affine flow, absence of flow induced ordering) inherent in the theory. Also, a qualitative change 
of the flow curves enables us to bracket the glass transition temperature T c of the theory from 
above and (for the first time in simulations) from below. Furthermore, the structural relaxation 
time in the steady state behaves quite similar to the system viscosity at all studied shear rates and 
temperatures. 
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Soft glassy materials under shear exhibit a rich phe- 
nomenology. In the dilute regime, at temperatures cor- 
responding to the liquid state, forced Raylcigh scattering 
experiments 0] show an increase of diffusion constant 
upon shearing (shear thinning), distinct from Taylor dis- 
persion (displacement of particles in the flow direction 
(x) as they move in the direction of shear gradient. This 
would give rise to < Ax 2 (t) >~ t 3 \) Q- At higher densi- 
ties, experiments show evidence for shear thinning due to 
the presence of freely slipping two dimensional crystalline 
layers 0- 

On the other hand, studies of disordered suspensions of 
hard spheres show that shear thinning and shear melting 
phenomena may also occur in the absence of a crystalline 
structure Q . Similar observations have also been made in 
light scattering echo studies of (disordered) dense emul- 
sions @. Brownian dynamics simulations show that 
shear thinning in concentrated colloidal suspensions is 
related to the fact that, in the limit of low shear rates, 
the main contribution to the shear stress originates from 
the Brownian motion of colloidal particles and that this 
contribution decreases with shear rate 7 [|| . 

Recent theoretical progress, connecting nonlinear rhe- 
ology with glass formation, prompts us to study in detail 
by simulations the stationary states of a well established 
glass model under shear. We focus on the yielding 
behavior for constant shear rate of glassy states close to 
but also far below the glass transition temperature. 

Recently, Berthicr, Barrat and Kurchan studied nu- 
merically a driven spin glass and showed that shear thin- 
ning can be understood in terms of an acceleration of 
inherently slow system dynamics by the external drive 
@. Within this approach, the stress depends on the 
shear rate via a power law (no dynamical yield stress). 
The 'soft glassy rheology model' of Sollich et al. Q ex- 
tends the minimal "trap model" originally introduced by 
Bouchaud 01 m order to take into account the effect 



of an external drive. The theory contains a noise tem- 
perature, x, which controls the distance from the glass 
transition at x = 1. For 1 < x < 2, a power law decrease 
of the stress with applied shear rate is found, whereas in 
the jammed state (x < 1), a continuous onset of a dy- 
namic yield stress is p redicted, a + = 17(7 — ► 0) = 1 — x. 
Fuchs and Cates [llj, on the other hand, started from 
the well studied mode coupling theory (MCT) of the 
glass transition [T^j. They started from the idealized 
picture that, in a supercooled liquid, nearest neighbors 
of a particle form a cage which progressively solidifies 
eventually leading to a complete arrest of all particles as 
the glass transition is reached. The effect of shear then 
enters by the advection of density fluctuations. Fluctua- 
tions of a given length scale are advected towards progres- 
sively shorter length scales so that particles need explore 
smaller regions in order for density correlations to decay. 
The interesting prediction of a yield stress discontinuity 
at the ideal glass transition was made. A related MCT 
approach to the fluctuations around the steady state has 
recently been proposed by Miyazaki and Reichman |l3j |. 
The issue of yield stress discontinuity, however, could not 
be addressed in that approach. 

Berthier and Barrat Q performed molecular dynamics 
simulation studies of the present model under a homo- 
geneous shear showing e.g. that, in a range of low shear 
rates, time-shear superposition and space time factoriza- 
tion theorems hold thus suggesting that generic prop- 
erties related to the glass transition "generalize" to the 
non-equilibrium situation of a homogeneous shear. How- 
ever, due to a rather limited range of shear rates, results 
presented in Ref . Q did not allow a clear answer whether 
the present model exhibits a yield stress or not. 

In this paper, we focus exactly on this aspect, namely 
an analysis of the dynamic yield stress and its behavior at 
the ideal glass transition. For this purpose, we performed 
large scale molecular dynamics simulations of a generic 
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glass forming system first introduced by Kob and Ander- 
sen 0- The model consists of a 80:20 binary mixture of 
Lcnnard-Jones particles (whose types we call A and B) 
at a constant total density. A and B particles interact 
via a Lennard- Jones potential. Standard (dimensionless) 
parameters will be used as in Refs. 0, El • 

Ten independent samples, equilibrated at a tempera- 
ture of T = 0.45, serve as starting configurations for all 
simulations reported here. The temperature is set to the 
desired value at the beginning of shear motion. The shear 
stress is computed using particle positions and velocities 
via the Irving-Kirkwood formula |15| . Only samples cor- 
responding to the steady state, i.e. to strains larger than 
100%, are taken into account (previous studies of the 
stress-strain relation of the same model showed that the 
initial transient behavior is limited to strains below 50% 
[T3j). Depending on the desired accuracy of the stress- 
data, the length of simulations was varied between 380% 
and 760% strain. Equations of motion are integrated us- 
ing a discrete time step of dt = 0.005. 

Recently, it was found that the present model may ex- 
hibit shear-localization in the glassy state if the shear 
rate is imposed by usin g a conventional Couette cell with 
moving atomistic walls In the present analysis, how- 
ever, we are interested in effects of a spatially constant 
shear rate, a basic ingredient of all theories briefly ad- 
dressed above. Therefore, we do not use atomistic walls 
but apply the so called SLLOD algorithm combined with 
Lees-Edwards boundary condition ^f|. With this sim- 
ulation method, we do indeed observe a linear velocity 
profile in all studied cases (not shown). 

As shear may, at least in principle, change the static 
structure of the system, we examine this by computing 
pair distribution functions along x (flow), y and z (shear 
gradient) directions separately. Figure Q illustrates re- 
sults on g(x), g(y) and g(z) (A-A correlations) for three 
characteristic temperatures T = 0.2 (deep in the glass), 
T = 0.42 (close to T c ) and T = 0.6 (supercooled state) 
indicating the absence of long range order in all stud- 
ied cases. Surprisingly, even though the stress changes 
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FIG. 1: A-A pair correlation functions along x (flow), y and 
z (shear gradient) directions, for three characteristic temper- 
atures: T = 0.2 (glass), T = 0.42 (close to T c ) and T = 0.6 
(supercooled state). The shear rate is j — 10~ 4 . For clarity, 
data at T = 0.42 (T = 0.2) are shifted upwards by 1 (2). 



by more than a decade at this 7, the local structure 
varies little with shear, remains amorphous, and (almost) 
isotropic. Similar observations are also made using the 
B-B and A-B pair correlations. 

Figure |21 shows simulated steady state shear stresses as 
functions of shear rate (viz. 'flow curves') for tempera- 
tures ranging from far above to far below the glass transi- 
tion temperature of the model. Note that shearing allows 
us to access stationary states at very low temperatures. 
As apparent from the change of the curvature of the flow 
curves (S-shaped without extended horizontal piece at 
high, horizontal piece merging into upward curvature at 
low temperatures) the system response changes qualita- 
tively around a temperature T c , which MCT identifies 
as (ideal) glass transition temperature. For T > T c the 
stress becomes proportional to shear rate as 7 approaches 
zero (linear response). At temperatures below T c , how- 
ever, a stress plateau forms in the low 7-regime. The 
observed qualitative agreement with the MCT scenario 
prompts us to test its predictions in more detail. Impor- 
tantly, the qualitative change manifest in the u[fi) curves 
enables us to give upper and lower limits for T c without 
any theoretical analysis: we conclude 0.34 < T c < 0.45. 
Fits with the schematic F^-model of MCT support this 
estimate, and achieve to describe the flow curves for all 
studied temperatures by adjusting two global parame- 
ters (a scale for 7 and one for a), and the parameter 
e at each temperature measuring the distance to the 
ideal glass transition 0|. The fit gives T c = 0.4 and 
°c~th = 0-19. We thus conclude that the basic rheologi- 
cal features of our model are well described within simple 
schematic models in the framework of the idealized MCT. 
Relaxation channels not contained in the idealized MCT 
(so-called 'hopping effects' 0, which may be the origin 
for the deviations from theory at very low 7 in Fig. |2J) 
can not falsify our bound for T c , because the qualitatively 
different shapes of the flow curves are the characteristics 
of the fluid or glassy states within MCT. 

The stress plateau is best developed for tempera- 
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FIG. 2: Simulated shear stress (symbols) compared to theo- 
retical predictions (solid lines) for various temperatures (T = 
0.01 to T = 0.60 from top to bottom). The inset shows the 
fitted separation parameter, e, versus T. 
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FIG. 3: Determination of dynamic yield stress and its tem- 
perature dependence (see text). 

tures deep in the glassy phase extending over about two 
decades in shear rate. Its onset is shifted toward pro- 
gressively lower 7 as the temperature is increased toward 
T c . This makes an estimate of the dynamic yield stress, 
cr + (T) = cr(T;7 — ► 0), a difficult task for temperatures 
below but close to T c . Nevertheless, an estimate of a + (T) 
is interesting because it highlights the anomalous weak- 
ening of the glass when heating to T c . Testing the MCT 
predictions below T c has previously not been possible in 
simulations because of problems to reach the equilibrated 
or steady state at sufficiently low shear rates. Our esti- 
mate is obtained by comparing the steady state shear 
stress for the two lowest simulated shear rates, namely 
7 = 1(T 5 and 7 = 3 x 1CT 6 . As shown in Fig. at 
temperatures below T = 0.38, practically the same shear 
stress is obtained for both choices of 7 indicating the 
presence of a yield stress plateau. 

For T — T c , we make use of theoretical predictions 
based on the F^ 2 -model For not too low shear rates, 
the flow curve takes the form of a generalized Hershel- 
Bulkeley constitutive equation, 

a = v+ (1 + \i/%\ m + c 2 |7/7*| 2m + c 3 |7/7*| 3m ) • (1) 

Here, 7* is an upper limit where this expansion holds. As 
the simple F^-model is found to describe the rhcological 
properties of our system rather well (see the discussion 
of Fig. El), we set the parameters, c 2 = 0.896, C3 = 0.95 
and m = 0.143 as obtained from an analysis of the model 
[T^ . We then apply a fit to Eq. ^ with cr+ and 7* being 
the only fit parameters. This gives cr+ sim = 0.15±0.01 
and 7* = 0.0045 ± 0.0008 the latter being close to our es- 
timate of the window, where MCT can describe the flow 
curves. At higher shear rates, 7 > 7,, we expect micro- 
scopic effects to dominate the stress. As shown in Fig- EI 
while (7 + (T) weakly varies with T at low temperatures, 
it steeply drops as T approaches T c , signalling the glass 
transition. The yield stress follows well the MCT-square 
root law 113 , cr+(T) - a+ oc |1 - T/T c |°- 5 . 

Note that the critical temperature of T c = 0.4, used in 
order to obtain best agreement between the theory and 
simulations, is slightly lower than the estimate T c = 0.435 
obtained from the analysis of the equilibrium dynamics 
of the system 0- This discrepancy is possibly related 



to the hopping effects observed in the density correlation 
functions closely above T c A closer analysis of this 

aspect requires understanding of hopping effects under 
shear and is beyond the scope of the present report. 

Next we address the close connection between observed 
non-linearity in the flow curves and the structural relax- 
ation as reflected in the incoherent scattering function, 

?«(*) =< E 4 =iexpfe(^(i) -ife(0))] > /N A . Here, N A 
is the number of A-particles, q y is the y-component of 
the wave vector and j/i is the y-component of the posi- 
tion of i-th particle. We set q y = 7.1 to the inverse of the 
average intcrparticle distance (or, equivalcntly, the max- 
imum position of the static structure factor) . The above 
definition of 3> 9 (i) takes into account the fact that the 
imposed flow is in the x-dircction and the shear gradient 
in the z-direction, and eliminates the advection of parti- 
cles with the flow. Note that the above choice of A type 
is arbitrary. Similar results are also obtained for type B 
particles (not shown). 

Figure H illustrates $ 9 (t) for T = 0.3 (glass) for all 
studied shear rates. First, & q (i) exhibits the two step 
relaxation process typical of supercooled liquids: A short 
time decay to a plateau, f q (characterizing the 'solidity' 
of the system at length scale 1/q 01) followed by a final 
decay to zero at much larger times. The long time decay 
of &q(t), on the other hand, is clearly dominated by the 
imposed shear. Its shape does not depend on the shear 
rate which only sets the time scale, T q ~ I/7. Note that 
Fig. 01 and the stationarity of the correlators, verifies 
that we achieved ergodic stationary states as expected 
from theory. 

Figure [5] depicts a = T q j (this definition of a stress 
goes back to Maxwell) for q = 7.1 versus 7 for various 
temperatures. Here, r q = J °° dt$> q (t) is defined as av- 
erage relaxation time. Where r q ~ I/7 holds, a plateau 
in a follows. Note the striking similarity between the 7- 
dependence of the shear stress (er; Fig. |2J) and that of a. 
As suggested by the inset of Fig. in the whole range 
of studied shear rates and temperatures, a /a changes at 
most by a factor of two whereas the shear stress varies 
by two orders of magnitude. Thus, both the 7 and T- 
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FIG. 4: Steady state incoherent scattering function, $ 9 (f), of 
a glass (T = 0.3) measured at various shear rates as indicated, 
and for q = 7.1. The quiescent $9(7 = 0) after a waiting time 
of t w = 10 is also shown. 
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FIG. 5: a = r g 7 versus shear rate for various temperatures 
ranging from the glassy phase to the supercooled state. The 
inset shows the ratio of this (Maxwell type) stress to the real 
shear stress for T = 0.2, 0.3, 0.4, 0.45, 0.525 and 0.6. 

dependence of steady state shear stress is mainly deter- 
mined by that of a = jr q . This observation is quite 
significant as it supports theoretical approaches where, 
even beyond the Newtonian-regime.the shear viscosity 
is simply taken as a relaxation time 0, IT^ | , or where this 
relation holds as approximation 

In summary, large scale molecular dynamics simula- 
tions have been performed in order to investigate the ex- 
istence and temperature dependence of the dynamic yield 
stress, cr + , for a 80:20 binary Lcnnard- Jones model glass 
first proposed by Kob and Andersen 0- Our data do 
indeed support the existence of a dynamic yield stress in 



the glassy phase as underlined by stress plateaus extend- 
ing over about two decades in shear rate. Let us men- 
tion recent experiments on the rheology of dense colloidal 
dispersions [18J which also find finite a + . Furthermore, 
the temperature dependence of o~ + follows the predicted 
anomalous weakening close to the glass transition with a 
finite critical yield stress of cr+ = 0.17 ± 0.02 as a com- 
promise between <r^ th = 0.19 and cr^ sim = 0.15. The flow 
curves allow for the first bracketing in simulations of the 
critical temperature of MCT from below. Irrespective of 
hopping effects neglected in the employed MCT, we can 
conclude T c > 0.34. 

Furthermore, a generalized stress (er; Fig. |5J) as the 
product of the shear rate and the structural relaxation 
time was determined for all temperatures and shear rates. 
a exhibits exactly the same qualitative features as the 
real shear stress thus emphasizing the close connection 
between the structural relaxation and the rheological re- 
sponse. Noting that a priori no relation a to a exists 1 1 1| . 
this finding corroborates theoretical approaches which 
make use of the relaxation time as a sort of viscosity. 
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